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We put forward a simple procedure for extracting dynamical information from Monte Carlo simula- 
tions, by appropriate matching of the short-time diffusion tensor with its infinite-dilution limit coun- 
terpart, which is supposed to be known. This approach -discarding hydrodynamics interactions- 
first allows us to improve the efficiency of previous Dynamic Monte Carlo algorithms for spherical 
Brownian particles. In a second step, we address the case of anisotropic colloids with orientational 

fSJ degrees of freedom. As an illustration, we present a detailed study of the dynamics of thin platelets, 

T-H with emphasis on long-time diffusion and orientational correlations. 

o 

C\j PACS numbers: 

OQ I. INTRODUCTION 

(~| Monte Carlo simulations (MC) provide a powerful method for calculating the thermodynamical averages of phys- 

O ical quantities of many-body systems and have been employed to study the equilibrium properties and phases of a 
2^ large variety of physical systems. Such approaches rely on the generation of a Markov chain -a stochastic set of 

configurations- with appropriate sampling of phase space. In the most common variants of MC, such as Metropolis 
algorithms, trial moves are accepted with a certain probability that satisfies detailed balance with respect to the de- 
ed sired Boltzmann distribution [U [2] . This method is frequently used to obtain average values of macroscopic quantities 
c/) in equilibrium thermodynamic ensembles: MC methods were initially devised to study static properties. In some case 
^_i though, where mesoscopic degrees of freedom interact with microscopic ones, the ensuing separation of time scales 
2 allows to replace the microscopic details by a noise, in the spirit of a Langevin equation. A suitably chosen MC scheme 

^ -stochastic in nature- then allows to study dynamical features of the mesoscopic degrees of freedom, meaning that 
JL( correlations between successive configurations in the Markov chain of the MC simulation can be interpreted in terms 

of the dynamic correlation functions. The legitimacy of this approach stems from the coincidence of the Fokker-Planck 
Q equations of the original system governed by the aforementioned Langevin-type dynamics, and of the fictitious MC 
Q dynamics [SHSj- 

'"^ Recently, there has been a rise of interest to employ MC simulations to study the dynamics of colloidal suspensions 

. [HinilT]. On the time-scale that momenta correlations have decayed, colloids undergo diffusive motion as a result of 

►>. collisions with solvent molecules. Therefore, at these time-scales the stochastic dynamics generated by MC algorithm 

^^ seems to be more appropriate compared to the deterministic Newtonian dynamics where solvent is omitted. Dynamic 

CK (sometimes known as Brownian) Monte Carlo (DMC) algorithms in which only single particle moves with sufficiently 

On small displacements are allowed, reproduce the real dynamics for times larger than the time-scale of momenta relax- 

\l^ ation. In such a case, DMC becomes equivalent to Brownian dynamics (BD) simulations [5HT^: a BD algorithm is 

ff^ also stochastic, with integrated out momenta and positions evolving with overdamped Langevin dynamics [1]. The 

<^) advantage of DMC over BD is that it is easily adaptable to systems with non-differentiable (hard) potentials. Al- 

CN though an "event-driven" variant of BD technique has been developed [13] to deal with such types of interactions, it 

\ I is computationally more cumbersome and expensive. Hence, studying the dynamics of hard particles with the DMC 

^ scheme seems to be an efficient route, provided that an accurate mapping between the Monte Carlo time-step and 

k>l ^^^ physical time is worked out. Achieving this goal is the main motivation of the present work. 

rN The significance of establishing the matching of time-scales is justified by the recent increasing use of DMC for 

j^ studying dynamics of various systems |14H17) . Recently, it has been proposed that equating the square of amplitude 

of MC displacement scaled with acceptance probability with infinite-dilution limit diffusion coefficient provides a 

good estimate of the physical time for spherical particles [B]. The applicability of this proposition is justified through 

agreement of BD simulations results with those of DMC [51[7]- Furthermore, these studies show that scaling the Monte 

Carlo time step with acceptance probability allows one to extend the limit of validity of DMC to relatively larger 

displacement amplitudes corresponding to acceptance probabilities -the fraction of accepted MC attempted moves- 

significantly smaller than 1. Here, we propose an alternative physically motivated approach for mapping the MC time 

to physical time. It allows us to push the limit of applicability of DMC to even larger displacement amplitudes (smaller 

acceptance probabilities). Our scheme is based on equating the short-time self-diffusion extracted from simulations 

directly with the infinite-dilution diffusion coefficient. Particular attention is paid to anisotropic particles. At variance 



^ 

a 



with previous approaches that did not consider the anisotropy of the short time diffusion tensor [7] , we have taken 
into account the important couphng between orientational and translational degrees of freedom. In all what follows, 
the various short-time diffusion constants are supposed to be known, and implicitly account for the presence of an 
underlying solvent. 

The rest of the paper is organized as follows. In sectionlll] the method is presented, for both spherical and anisotropic 
particles. In particular, we discuss the relation between the amplitudes of translational and rotational moves, essential 
to achieve a physically consistent diffusive process. In section |III[ we study the convergence and self-consistency of 
DMC simulations as a function of displacement amplitude for both spherical and disk-shaped particles, and compare 



our approach with previous investigations. As an illustration, the method is employed in section |IV| to explore the 
development of long-time translational diffusion of infinitely-thin hard disks (platelets) as a function of density. We 
find that upon increasing the density deep in the nematic phase, the long-time diffusion becomes anisotropic, and that 
in contrast to an impeded diffusion in the nematic direction, the transverse diffusion of disks is enhanced. Concluding 
remarks close the paper with section \V\ 

II. METHODOLOGY 

We start by describing the dynamic Monte Carlo algorithms used for both spherical and anisotropic particles. In 
each case, we discuss the procedure for matching of the time scales. We then introduce the model systems and provide 
the simulation details. 

A. DMC algorithm for spherical particles 

Colloids suspended in a solvent undergo overdamped Brownian motion with diffusive behavior, for large enough 
times compared to the momentum relaxation time t|^. The latter quantity is set by the colloids mass AI and the 
translational friction coefficient 74 that depends on the particle size and shape, its value being Sttt^ct for spherical 
objects of diameter a with stick boundary conditions: r|^ = M/^t. The resulting mean-squared displacement (MSD) 
of non-interacting colloids varies linearly with time for t ^ t\^j, with a slope given by the infinite-dilution diffusion 
coefficient Dq = kBT/^t- However, for interacting colloids in non-dilute suspensions, different diffusion processes 
should be distinguished, namely short-time D)^ and long-time D\^ diffusion. The distinction requires the introduction 
of the Brownian time-scale tb defined as the time required for an isolated colloid to diffuse over its diameter a^ i.e., 
tb = ct^/(6Dq). For relatively short times, larger than t|j but smaller than the Brownian time-scale, the colloids 
influence each others motions indirectly through the solvent flow field in which they move. These solvent mediated 
hydrodynamic interactions may affect the short-time diffusion. If one ignores the hydrodynamic interactions, as in 
the subsequent analysis, the short-time diffusion is that of infinite-dilution diffusion coefficient D\ = Dq [THlIin]. On 
the contrary, the long-time diffusion Z?^ that is defined for i ^ tb is mainly determined by the direct interactions 
between colloids [13120] ■ For typical colloids with diameters in the range 10 nm-1 /im, we have tb/t\,j :^ 1 with well 
separated diffusive regimes. 

Now, consider a MC procedure discarding hydrodynamic interactions, where each of N interacting spherical particles 
in the simulation box is shifted by a random displacement chosen in the interval [—SI, 61] along each Cartesian 
coordinate. The moves are accepted according to the Metropolis algorithm [2]. Such a simulation mimics the Brownian 
motion of the colloids for time-scales t ^ t'Ij. One expects that the mean-squared displacement (MSD) of a particle 

after n cycles (Ar^(n)) — l/A^ X]i=i(l'^i('^) ~ ^i(O)P) varies linearly with the number n of MC cycles, for both small 
values of n corresponding to the short time regime, and large n, albeit with a different slope. The MSD for sufficiently 
small n is governed by the infinite-dilution diffusion coefficient Dq, i.e., (Ar^(n)) w 6DQ{n5t), provided that the 
amplitude of MC move SI is chosen sufficiently small, i.e., SI <^ a where St is the physical time interval that each 
MC cycle corresponds to. As a result, we impose that the relation between the MC clock and the real time St can be 
obtained from the slope of MSD in the small n limit, i.e. the short-time diffusion of MC simulation. 

SI ^ ^^ (A^ ^^^ 

Tb "->o na'^ 

where lim„_5.o with ji an integer refers to the limiting behavior of the MSD slope at small n. It has been noted that 
for sufficiently small SI, (Ar^(l)) = ASP where A is the acceptance probability of the MC scheme [6'. It was therefore 
suggested that the time-scale corresponding to a MC cycle can be obtained as: St = ASP /6Dq . Scaled in terms of 
Brownian time, this equation can be written as 

^ = a'4 (2) 

Tb cr-' 



which provides an alternative route against which our approach will be tested in section III This ^-rescaling procedure 
will lead below to the variant denoted Va, while the diffusion-matching will be referred to as variant Vd- We will 
show that, quite expectedly, both methods become equivalent, for sufficiently small SI. However, enforcing (IT]) instead 
of pi) allows us to extend the limit of applicability of DMC towards larger values of SI. 



B. DMC algorithm for anisotropic particles 

For anisotropic particles with orientational degrees of freedom and depending on the shape, the diffusion in some 
directions is favored over some others, leading to the coupling of translational and rotational motions [211. Henceforth, 
for a meaningful dynamics, due account should be taken of the anisotropy of diffusion in the body-frame. Similar 
to the translation-only case, simulations based on DMC can produce the correct dynamics of the rotational degree 
of freedom for time scales larger than the damping time of angular velocity tI,j = Ir/"fr where Ir is the moment of 
inertia and 7^ is the rotational friction coefficient. The exact value of t^,j depends on the size, shape of the particle 
and the axis of rotation under consideration, but its order of magnitude is the same as r^. On the other hand, the 
time-scale for orientational relaxation t^ is given by the inverse of the infinite-dilution rotational diffusion coefficient 
Dq = ksT/'^r leading to Tt/tJ^j = ksT/Ir 3> 1 for sufficiently large colloids. 





FIG. 1: A schematic drawing presenting a) an example of an axially symmetric object with its orientation vector Ui together 
with its body frame defined by Ui and two mutually perpendicular vectors rixi and riyi. b) rotation of orientation vector with 
an angle 56. 



To illustrate the DMC implementation for anisotropic colloids, we focus here on axially symmetric particles, although 
the generalization of our method to less symmetric objects is straightforward. We consider a system of N thin 
disks. Each of them can be identified by its center of mass position r^ and its orientation unit vector u^, that is 
taken along the symmetry axis. The translational (resp. rotational) diffusion tensor in the body frame is diagonal 
and consists of one coefScient Dqm (resp. DL,) for the direction parallel to Ui and two identical coefficients Z?qj^ 
(resp. Dqj_) for perpendicular directions. In the following, we describe two variants of the same MC algorithm, 
that consists of simultaneous translational and rotational displacements enforcing proper symmetry of the diffusion 
tensors. We ignore the rotations around the symmetry axis of the particle, as they cannot be detected in most of 
experiments; only rotations around axes perpendicular to symmetry axis are considered, i.e rotations of the orientation 
vector, characterized by the same coefficient DJ^j_ which for simplicity we denote as Dq. The first variant, Va, bears 
similarities with the approach presented in [7] but differs in two respects: first, no fine tuning of translational and 
orientational moves are required here; second, we fully resolve the anisotropic dynamics at short times, which is an 
essential prerequisite for the study of phases with orientational order, such as nematic phases. The second variant is 
reminiscent of the method put forward above for spherical particles. 

We start by considering the limit of small MC increments, for which the move acceptance probability is close to 
unity. A general attempted translational move can be described by Sfi = Sxrixi + Syriyi + Szui, where na;^ and fiyi 
are two mutually perpendicular arbitrary unit vectors that are normal to Ui as well, Sz is a random number in the 
interval [—SI, SI] and Sx, Sy are random numbers in the interval [—SI' , SI']. For sufficiently small SI and SI' , the average 
MSD in the direction parallel to Ui is {Sz'^) ~ SP/3 and along perpendicular directions is (fe^) = {Sy'^) — Sl'^/3. 



Imposing the symmetry of the diffusion tensor provides us with a relation between SI' and SI , i .e. 

SV 




SI M/Z?*|| 



(3) 



In addition, in the short-time diffusion regime, we should have {Sz'^) — 2DQut, so that the time increment corresponding 
to one single MC cycle reads 

(fe^(l)) 
2Di 



^t = '^^^^ (4) 



'oil 

If the orientation distribution for the ensemble of disks is isotropic, a restrictive assumption, we have {Sr^) = 
2Do±Sl'^ + Dq^^SP. Defining the average translational diffusion as Dq — {2Dqj_ + DL.)/3, the MSD can then be 

simphfied to (Sr^) ^ ^Sf. 

The change of orientation can be seen as a random rotational displacement with an angle SO in the interval [0, 5a], 
such that u^'^'^.Ui = cosS6, see Fig. Il|for an illustration. Such a rotation can be achieved as follows [2 [52]: we 
generate a unit vector uj with an isotropic random orientation and obtain the new orientation vector as 

^ne^ = AA[(1 - Sa)u, + Sau'i] (5) 

where M ensures proper normalization. From this, one can calculate the correlation between the new and initial 
orientation vector in terms of Sa: 

{uT'"-Ui) = F{Sa) = g/ ^ ^ ^ ^2 -5" < 0-5 (6) 

~ 1 - Sa^/i + 2/3fo^ Sa < 1 

where the average is taken over all the possible orientations of the random vector u[ . In addition, for the physical 
Brownian system under study, we have (u(t).w(O)) — exp(— 2D5i), |23] so that we get from equation (|6| that 

ln(F(M) (,) 

In the limit of small Sa, the mean-squared angular displacement for diffusion of orientational vector can be obtained 
from Eq. ^■. \{u{Sa) - u(0))p = {SO"^) ~ 2Sa'^ /i (for Sa < 1) that should be equal to ADlASt. The last step is to 
enforce consistency of time scales, equating the two relations for St, i.e. Eq. Q and Eq. (l7|. In doing so, we obtain 



^t-^r^ (8) 




'oil 
and the following constraint between the amplitudes of translational and rotational moves: 

Sl = ^~2,HF{Sa))J^ (9) 

V oil 
In the limit of smaller Sa, this simplifies into: 

x„. I rtr 

(10) 

We should now take due account of the fact that A ^ 1. In the first variant of the approach, Va, it is assumed that 
the physical time increment is slowed down by rejected moves, so that the generalized Eqs. (It]) and (IS]) read: 

SP 
St^A-^ (11) 

fin* 

^^^_^ln(F(M) (12) 

21)5 

Therefore, once the amplitude of orientational moves Sa has been chosen, the amplitudes of translational moves SI 
and SI' follow from Eq. ^ and Eq. (l3|. Physical time is given by (11) or cquivalently (12 1, where the acceptance 



probability A is computed on the fly in the simulation. Alternatively, for variant Vd, we again impose ([3]) and (|9|, 
but determine the physical time scale by imposing that ((5z^(l)) computed in the simulation, coincides with 2DL,t . 
Before illustrating the applicability of our DMC algorithm, we provide in the following section some details con- 
cerning the systems simulated. 



C. Model systems and simulation details 

In the remainder, we investigate the dynamics of two model systems by means of DMC simulations: hard sphere 
colloids (system A) and infinitely thin disks (system B) with diameters a. We take into account the direct hard-core 
interactions by choosing at random a particle and generating a random trial MC move (including the rotational move 
for disks). We then reject the displacements that lead to an overlap with neighbors ^2] . 

The first system consists of N particles in a cubic box of length L, with periodic boundary conditions. We took 
N = 1024 and the simulations were performed for volume fractions $ = ttct^ / {6N L^) in the range 0.05 — 0.5. The 
starting configuration was that of a BCC crystal melted by an equilibration run of 2 x 10^ MC cycles (one trial move 
per particle). The production runs for calculating the mean-squared displacements consisted of 1 — 5 x 10^ cycles, 
depending on the volume fraction and the amplitude of translational displacements 0.01 < Sl/a < 0.1. 

The second system consists of 500 disks again in a cubic simulation box with periodic boundary conditions. For each 
reduced density p* = Na^ /L^, a first equilibration run of 5— 10 x 10'' MC cycles starting with an initial configuration 
of disks on an FCC crystal with parallel orientations. The production runs for calculating the translational mean- 
squared displacements and orientational correlations were in the range of 10^ to 10^ cycles, depending on the density. 
The displacements amplitudes were in the range 2 x 10"'' < Sl/a < 0.2. and 0.005 < 6a < 0.324. The infinite-dilution 
translational and rotational diffusion coefficients of disks used in the simulations are 

^011 = ^ and D„^ = ^^ (13) 

giving an average diffusion coefficient of Dq = ^^. On the other hand, we have for the rotational diffusion 

^0. = !^. (14) 

These results are obtained from the general formula of diffusion coefficients of oblate spheroids (Ml [21] in the limit of 
vanishing length of semi-minor axis. Having described the methodology and simulation details, we present below the 
results of our DMC simulations. 

III. ASSESSMENT OF THE DYNAMIC MONTE CARLO SCHEME 

In this section, we present our DMC results for hard sphere self-diffusion, and then turn to thin disks. We compare 
and discuss the two different procedures for mapping MC time, i.e. rescaling with acceptance probability or directly 
matching short-time dynamics with infinite-dilution diffusion tensor. 

A. Dynamics of spherical colloids 

We start by discussing the time-scale matching in a colloidal suspension of hard spheres. Figure [2] shows the time 
dependence of the mean square displacement at a relatively high volume fraction $ = 0.5, for different values of MC 
displacement amplitude S = 51/ a. In Fig. 2k, the physical time is obtained from (5^ scaled with acceptance probability, 
i.e., I/tb — nAS^ as suggested in Ref. [S . As reported in 6J, such a procedure leads to a decent data collapse, the 
goal being to obtain results that do not depend on 6. In this respect, the collapse is only partial, see e.g. the (5 = 0.1 
data that do not completely superimpose to those for 6 = 0.01. 

In Fig. I2b, we have obtained the physical time from the alternative method leading to Eq. (fTl). The graph shows 
that equating the short-time diffusion from MC with Dq directly, allows a better collapse of MSDs for larger values 
of (5 (where acceptance probabilities are smaller). We emphasize that for the small 6 = 0.01, the two approaches yield 
the same results. However, with the present proposal, we can employ the DMC algorithm for relatively larger values 
of the increment S. 

As discussed in the methodology section, the rationale behind our method is the fact that at short-times (t <C t^), 
the particles diffuse freely with diffusion coefficient Dg, while for long enough times {t ^ tb), the MSD crosses 
over from free diffusion to a slowed-down motion characterized by long-time diffusion _D^. This is illustrated in Fig. 
[Sk, where the reduction of diffusion coefficients is ten-fold. To further test the reliability of our method, we have 
plotted in Fig. [3)3 the long-time diffusion coefficient, extracted from the slope of the MSD curve at long times, as a 
function of volume fraction $ for different sampling amplitudes. For comparison, we have also included the results 
obtained from scaling with A for 6 = 0.1 and Brownian dynamics results taken from reference [26j . The simulation 
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FIG. 2: Mean-squared displacement obtained from DMC simulations for a hard sphere system of volume fraction $ = 0.5, 
as a function of time. The curves correspond to three different values of displacement amplitude 5 = ilja — 0.01,0.07 and 
0.1. The corresponding acceptance probabilities A are: 0.84, 0.317 and 0.202, respectively, a) t/rs is obtained by scaling with 
acceptance probability A using Eq. ([2|. The lower dashed curve shows the MSD data for S = 0.1 without ^-rescaling. b) I/tb 
is obtained from Eq. al\ by matching the MC short-time diffusion coefficient to Dq. For comparison, the curve with S — 0.01 
and yl-rescaling [i.e. that shown in panel a)] is also plotted. 



data from event-driven BD |TB] are fully consistent with our DMC results (not shown). While for 6 = 0.01 the two 
approaches are equivalent (as illustrated in Fig. [2t)), some discrepancy is visible for 6 — 0.1 and <& > 0.35. We 
conclude here that the yl-rescaling fares somewhat worse. For completeness, we also have displayed the theoretical 
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FIG. 3: a) The mean square displacement obtained for <1> = 0.5 and 5 = 0.002 clearly demonstrates a short-time slope of 6Dq 
(unity in scaled units) while at long time, the MSD grows with a reduced slope given by D^. b) Long-time diffusion coefficient 
of hard sphere colloids without hydrodynamic interactions, obtained from DMC simulations with two methods, matching of 
short-time diffusion from MC with Dq (solid squares 5 = 0.01 and triangle, S — 0.1) and with A-rescaling (open triangles, 
S = 0.1). Also shown are hard sphere BD results from reference [26] (open circles). The short-dashed curve is for the Tokuyama 
and Oppenheim formula 27 , that is used in reference [261 to obtain the ratio of long to short time diffusion. 



results of Tokuyama and Oppenheim for the ratio of long-time to short-time diffusion D\^/ Dg obtained for spherical 
particles when hydrodynamic interactions are accounted for. Good agreement is found with the present simulation 
data that discard such interactions, consistently with the findings of Ref. [5D]. 



B. Dynamics of thin colloidal disks 
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FIG. 4: DMC simulations results of thin platelets at a low density of p* — 0.1, obtained for a displacement amplitude of 
S = 5l/a = 0.005, for which the acceptance probability is 0.9988. a) The MSD as a function of I/tb = nSl^ agrees well with 
the line of slope unity b)The first and second order time orientational correlations versus t/ri^ = 2nDgSa^ /3 show a good 
agreement with exp(— f/r/''^). 

We now turn to the DMC simulations of thin disks. First, we investigate the time behavior of MSD and orientational 
correlations for a low density system, and compare these results with the theoretical expectations for a freely diffusive 
particle. As can be seen in Fig. |4^, the particles diffuse, as they should, with the same diffusion coefhcient at short and 
long times. To quantify orientational dynamics, it is customary to define the correlation functions {Pi{ui(t) ■ Ui(0))), 
where Pi is the l—th order Legendre polynomial. For a colloid in diluted conditions, these orientational time correlation 
functions decay exponentially with a relaxation time rj^^ — 1/{DqI{1 + 1)) [23, 25 . Of particular interest among the 
correlation functions are those associated to Pi and P2, that are related to the dielectric properties of polar liquids 
and to the scattering of depolarized light, respectively [3S]. In panel (b) of Fig. [4J we have plotted both the first and 
the second order orientational correlation functions versus time. These functions show a very good agreement with 
their analytical infinite dilution counterparts. It should be noted that with the parameters chosen in Fig. HI where 
the acceptance ratio is close to unity, variants Va and Vd coincide. 

The next step is to explore the self-consistency of our two variants, where the time behaviour generated should 
be independent of the auxiliary parameters chosen for MC sampling. We have one such parameter, say da, from 
which the other relevant increments SI and SI' follow, see Eqs. (|9| and (Is]). We first analyze the behavior of the 
translational self-diffusion. In figure [SJ we have plotted the MSD of disks at a density p* — 2, that is below the 
density of the isotropic-nematic transition p*j^ = 4 [55] , as a function of physical time I/tb , for different MC sampling 
amplitudes. As discussed in section lllj we perform simultaneous translational and rotational moves. It appears that 
both variants Va [with results shown in panel a)] and Vd [results in panel b)] are satisfactorily self-consistent, with 
a proper collapse of data. Relatively large values of the sampling parameter are therefore acceptable, and provide 
results of a comparable accuracy as more demanding simulation with finer resolution. Upon closer inspection, it can 
be seen that variant Va shows a somewhat smaller dispersion of results that Vd ■ The analysis of orientational time 
correlations corroborates this conclusion, see Fig. |6] These conclusive tests allow us to study the density dependence 
of long time diffusion in a system of disks, and in particular the effect of a phase transition crossing. 

IV. ANOMALOUS DIFFUSION IN THE NEMATIC PHASE OF THIN DISKS 



Thin platelets undergo an isotropic-nematic transition upon increasing the density [22^, and it is interesting to see 
how the long-time translational diffusion and orientational relaxation are affected. To this end, we have performed 
systematic DMC simulations with sufficiently small displacement amplitudes and have obtained both long-time trans- 
lational self-diffusion coefficient and orientational relaxation time as a function of density. Fig. [7^ depicts the 
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FIG. 5: Mean-squared displacement obtained from DMC simulations for infinitely thin platelets of reduced density p* = 2, as a 
function of physical time. Different values of sampling amplitudes were used: 5 = ilja — 0.01, 0.05, 0.1 and 0.2 corresponding 
to Sa — 0.0225,0.109,0.196 and 0.324, respectively, and to the following acceptance probabilities : 0.93, 0.68 and 0.48. a) I/tb 
is given by variant Va b) Results of variant Vd ■ 
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FIG. 6: Same as Fig. |5] for orientational correlations. In both cases, the results are reasonably independent of the sampling 
parameters. 



long-time translation self-diffusion coefficient I?^. We find that increasing the density, D^ decreases up to the tran- 
sition point. However, in the nematic phase not only the short-time diffusion is anisotropic, but also, the long-time 
diffusion becomes anisotropic with respect to the nematic director. Upon further increasing the density, we observe 
that the diffusion coefficient in the direction perpendicular to nematic axis D*j^j_ grows while the parallel component 
D^ll decreases significantly. As can be seen from Fig. PZa, £'^_|_ approaches the free diffusion coefficient of disks 
Dqj_/Dq = 9/8 = 1.125 [see Eq. (13)] in the limit of very high densities. In Fig. ub^ we have plotted D^j^ versus 
nematic order parameter 5*. As demonstrated by this figure the more the disks become aligned, the larger is the 
perpendicular component of the diffusion in contrast to D^n that becomes very small: topological constraints due to 
the excluded volume constrain the disks to move in a caging slab of parallel neighboring particles. We also emphasize 
that for all results presented in this section, variants Va and Vd provide strictly identical results. 

At this point it is interesting to compare our DMC results with molecular dynamics (MD) simulations of thin disks 
[28j . where an anisotropic diffusion in the nematic phase has also been observed [21]. However, one should keep in 
mind that the model in this work is not equivalent to ours. In our DMC simulations, we mimic the presence of an 
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FIG. 7: a) Long-time average difFusion coefficient 73^ (in botli isotropic and nematic phases) and long-time diffusion coefficients 
in directions parallel D^n and perpendicular D^x to the nematic director for densities p* > 4. b) D^x ^-s a function of nematic 
order parameter. 



underlying solvent through the stochastic nature of MC moves. On the other hand, there is no solvent and hence the 
short-time difFusion is replaced by a ballistic regime in the MD approach of Ref . [28| . For a quantitative comparison, 



Fig. p^ shows the ratio -D^j^/D^,, as a function of density. It seems that the two models have a different limiting 
behavior at high densities, where the MD data exhibit enhanced anisotropy. 
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FIG. 8: Ratio D\^ID\^\ for hard disks as a function of density. DMC data are compared with their MD counterpart of Ref. 



We also have investigated the evolution of orientational time correlations with density, see Fig. [9^. We observe 
that the orientational time-correlations decay exponentially in the isotropic phase, while they become non-ergodic in 
the nematic phase and develop a plateau at long times whose value is equal to the square of nematic order parameter 
S. To quantify the development of relaxation time with density, we have fitted the orientational correlation functions 
in the isotropic phase with an exponential, and obtained the corresponding relaxation time for the first and the 
second-order correlations as depicted in Fig. [Sb. As expected, relaxation times grow with density upon approaching 
the isotropic-nematic transition. 
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FIG. 9: a) The second-order orientational correlation function for different densities of disks b) The first and second-order 
rotational relaxation times in the isotropic phase, as a function of density. 



V. SUMMARY AND CONCLUSIONS 

To summarize, we have presented a dynamic Monte-Carlo algorithm (DMC) for both spherical and anisotropic 
colloids. In each case, we have discussed the procedure for matching the Monte Carlo time scale with its physical 
counterpart. In the case of spherical particles, we found that matching the short-time diffusion from DMC to the 
infinite-dilution diffusion coefficient leads to a better convergence of results than acceptance rate based schemes, for 
relatively large values of displacement amplitude. A slightly better agreement of the long-time diffusion coefficient 
with Brownian Dynamics data available in the literature was thereby achieved. 

For anisotropic colloids, we presented two variants of the DMC algorithm that takes into account the anisotropy 
of short-time diffusion into account. As for spheres, one is acceptance-rate based (Va), and one relies on short- 
time diffusion matching (Vd). Both routes are new in their present formulation, although variant Va shares common 
features with the approach of Refs. [SI [7]. A key point is that the appropriate ratio of translational and rotational move 
amplitudes is enforced, which leads to the proper short-time diffusive behaviour. We have tested the self-consistency 
of both variants, that give similar results for a system of thin platelets in three dimensional space. The method 
was finally employed to investigate the evolution of the long-time diffusion coefficient and orientational correlation 
functions with density. The anisotropy of the long time translational diffusion tensor was characterized in the nematic 
phase. While diffusion along the nematic axis becomes small when nematic ordering is more pronounced, it is enhanced 
in the perpendicular direction. 
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List of symbols 

N: number of particles in the simulation box of size L; 

M: Mass of particles; 

(7 = 2R: Sphere or disk diameter; 

n: number of Monte Carlo cycles, where a cycle is defined as one MC move per particle; 

$ = 7rcr^/(6iV_L'^): volume fraction of spheres; 

p = N/L^: number density of thin disks; 

p* = pa^: dimensionless number density of thin disks; 

^tir)'- translational (rotational) friction coefficient; 

t\j — M l^t : time-scale for vi^hich momenta of Brownian particles have relaxed; 

t'^i — Ir/"fr'- damping time of angular velocity for rotational Brownian particles; 



Dq — kBT/"ft' infinite-dilution translational diffusion coefficient of spheres or average translational diffusion 

W3 + ^^l|/ 



coefhcient of anisotropic particles; for disks, we have Dq — 2Dqj_/3 + Di,J3. 



Dg] short-time translational diffusion coefficient; 

D\^: long-time translational diffusion coefficient; 

tb ^ (t^/(6-Dq): Brownian time-scale, required for diffusing over a distance equal to the particle size; 

Dqj_ = ksT/^^: infinite-dilution translational diffusion of an axially symmetric particle in the direction per- 
pendicular to the symmetry axis; 

£)*ll = ksT Iti- infinite-dilution translational diffusion of an axially symmetric particle in the direction parallel 
to the symmetry axis; 

i'^j^ : Long-time translational self-diffusion of axially symmetric particles in the direction perpendicular to the 
nematic axis; 

Dm : Long-time translational self-diffusion of axially symmetric particles in the direction parallel to the nematic 
axis; 

Dgj^ = Z?q: infinite-dilution rotational diffusion of axially symmetric particle in the direction perpendicular to 
the symmetry axis; t^ = 1/(2_D5): time-scale for relaxation of orientation vector; 
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„ll — h„^l^X. 



• Dgi, = ksT/yr — ksT/j:^: infinite-dilution rotational diffusion of axially symmetric particle in the direction 
parallel to the symmetry axis; 

• (Ar^): mean-square displacement (MSD); 

• ((5r^) = (Ar^(l)): mean-square displacement after one MC step; 

• (56"^): angular mean-square displacement; 

• A: acceptance probability; 

• SI: amplitude of the translational displacement 

• J = Sl/a 

• 5a: maximal amplitude of the rotational displacement 69; 

• dt: physical time interval corresponding to one MC cycle; 

• t[: the relaxation time of the orientational time correlation functions {Pi(u{t) ■ u(0))); 

• T^"^ ~ l/(^(^ + ^)Dq): the l-th order orientational relaxation time of an isolated particle; 

• S: nematic order parameter. 



